Kompresja

Pomocnicza funkcja rysująca obrazki

plotImg <- function(m) {
  par(mar=c(0,0,0,0), mai=c(0,0,0,0))
  dims <- dim(m)
  plot(c(0, dims[1]), c(0, dims[2]), type='n', xlab="", ylab="")
  rasterImage(m, 0, 0, dims[1], dims[2], interpolate=TRUE)
}

Otwieramy obrazek testowy

library(png)
# read img as matrix
img_png <- readPNG("./dog.png")
plotImg(img_png)

Zmieniamy format

Kompresujemy tylko warstwy RGB, kanał alpha pomijamy, ponieważ daje złe efekty. Nawet delikatne szumy na tym kanale istotnie pogarszają jakość całości. Wszystkie kanały układamy wzdłuż szerokości.

# drop alpha channel and align layers along width axis
imgwide <- matrix(img_png[,,1:3], nrow=nrow(img_png))
# save height of image
h <- nrow(imgwide)
plotImg(imgwide)

Obliczamy PCA

Wyłączamy centrowanie i skalowanie, ponieważ obrazki już są znormalizowane do [0, 1].

pca <- prcomp(imgwide, center=FALSE, scale=FALSE)

Tabela kompresji

Sprawdźmy ile razy mniejszy obrazek dostanejmy przy wypraniu pc wektorów. Rozmiar jest podany w ilości liczb.

compression <- data.frame(
  pc = 1:h,
  # size of compressed image is sum of sizes of two matrices required to recreate image: x, rotation
  size = sapply(1:h, function(k) length(pca$x[, 1:k]) + length(pca$rotation[,1:k]))
)
compression$size_ratio <- compression$size / length(imgwide)
knitr::kable(head(compression, n=20))
pc size size_ratio
1 4736 0.0017428
2 9472 0.0034857
3 14208 0.0052285
4 18944 0.0069713
5 23680 0.0087141
6 28416 0.0104570
7 33152 0.0121998
8 37888 0.0139426
9 42624 0.0156854
10 47360 0.0174283
11 52096 0.0191711
12 56832 0.0209139
13 61568 0.0226568
14 66304 0.0243996
15 71040 0.0261424
16 75776 0.0278852
17 80512 0.0296281
18 85248 0.0313709
19 89984 0.0331137
20 94720 0.0348565

Porównanie wyników

pcs <- c(1, 3, 5, 7, 10, 12, 15, 20, 25, 30, 40, 50, 75, 100, 150, 200)
# for each number of principal components we decompress image to test it quality
sapply(pcs, function(k) { 
  cat(paste0("Principal components: ", k, " Size ratio: ", compression[k, ]$size_ratio), "\n")
  img <- pca$x[,1:k] %*% t(pca$rotation[, 1:k])
  # normalise
  img[img > 1] <- 1
  img[img < 0] <- 0
  # split width into 3 layers
  dim(img) <- c(nrow(img), ncol(img)/3, 3)
  plotImg(img)
}) -> tmp

Principal components: 1 Size ratio: 0.00174282703030517 Principal components: 3 Size ratio: 0.00522848109091551 Principal components: 5 Size ratio: 0.00871413515152586 Principal components: 7 Size ratio: 0.0121997892121362 Principal components: 10 Size ratio: 0.0174282703030517 Principal components: 12 Size ratio: 0.0209139243636621 Principal components: 15 Size ratio: 0.0261424054545776 Principal components: 20 Size ratio: 0.0348565406061034 Principal components: 25 Size ratio: 0.0435706757576293 Principal components: 30 Size ratio: 0.0522848109091551 Principal components: 40 Size ratio: 0.0697130812122069 Principal components: 50 Size ratio: 0.0871413515152586 Principal components: 75 Size ratio: 0.130712027272888 Principal components: 100 Size ratio: 0.174282703030517 Principal components: 150 Size ratio: 0.261424054545776 Principal components: 200 Size ratio: 0.348565406061034

Sygnały

Ładujemy pomieszane sygnały

sgn <- read.csv("signals.tsv", sep="\t")

par(mfrow=c(2,2))
for (i in 1:4) {
  plot(sgn[, i+1], type="l")
}

Używamy FastICA do rozdzielenia sygnałów przy użyciu dwóch algorytmów.

Pierwszy algorytm wydaje się dawać lepsze rezultaty.

ica <- fastICA(sgn[2:5], n.comp=4, alg.typ="deflation")
par(mfrow=c(2,2))
for (i in 1:4) {
  plot(ica$S[, i], type="l")
}

ica <- fastICA(sgn[2:5], n.comp=4, alg.typ="parallel")
par(mfrow=c(2,2))
for (i in 1:4) {
  plot(ica$S[, i], type="l")
}